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INTRODUCTION 


The  idea  of  a semi-Maricov  process  was  proposed  in  1954  (references  1 and  2).  The 
semi -Markov  process  Is  similar  to  the  Maricov  process  in  that  both  processes  are  de- 
scribed by  a set  of  states  whose  transitions  are  governed  by  a transition  probability 
matrix.  The  semi-Markov  process,  however,  differs  from  the  Markov  process  in  that 
the  times  between  transitions  may  be  random  variables.  Further,  the  amount  of  time 
spent  in  any  state  after  entering  it  is  a random  variable  described  by  a probability  density 
that  can  be  a function  of  both  the  state  of  occupancy  and  the  states  to  v.rhich  transitions  can 
occur. 

The  statistical  time  behavior  of  the  semi -Markov  process  is  described  by  a set  of 
linear  integral  equations . The  solution  of  this  set  of  equations  yields  the  probabilities 
that  the  process  occupies  the  states  of  the  system  as  a function  of  time.  In  many  prac- 
tical cases  where  the  semi -Markov  process  is  complex  or  involves  many  states,  an 
analytic  solution  is  difficult  to  obtain  and  numerical  procedures  must  be  used. 

The  purposes  of  this  Research  Contribution  are  (1)  to  acquaint  the  analyst  with  the 
semi-Markov  process,  (2)  to  illustrate  that  use  of  the  semi-Markov  process  is  common  in 
systems  and  operations  analysis,  and  (3)  to  show  how  numerical  solution  of  the  associated 
equations  can  be  accomplished. 

The  general  formulation  of  the  semi -Markov  process  as  developed  in  references  3 and 
4 is  discussed  informally.  A more  formal  treatment,  however,  appears  in  references  5 
or  6.  Numerical  techniques  are  emphasized  since  algorithms  for  analyzing  semi -Markov 
processes  are  not  well  developed  or  well  known.  Algorithms  for  both  continuous -time  and 
discrete -time  processes  are  derived,  and  their  numerical  accuracies  are  discussed.  Two 
examples  are  worked  out  in  detail . 
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GENERAL  FORMULATION  OF  A SEMI-MARKOV  PROCESS 


BASIC  FEATURES 

Let  N be  the  number  of  states  of  a semi-Markov  process.  Let  h^j(t)dt  be  the 

probability  of  a transition  to  state  j between  t and  t+dt  given  that  the  last  transition 
was  to  state  i at  time  zero.  These  holding  time  functions  hy(.)  are  elements  of  the 

transition  matrix  H(.),  i.e.. 


H(t)  = (hj^(t)} 


1 « i.J  s N 


and  must  satisfy 


/dt  E h./t)  = 1,  1 « i « N . (1) 

0 j=l  J 

The  functions  h^j(*)  are  "dishonest"  density  functions  in  that  they  need  not  integrate  to 

unity,  although  they  are  otherwise  like  density  functions.  Equation  (1)  expresses  the 
Intuitive  requirement  that  there  is  unit  probability  that  the  system  will  be  in  one  of  the  N 
states  of  the  system  at  some  point  in  the  future,  given  start  in  i . 

When  the  system  makes  a transition  from  state  i to  state  J , the  transition  is  said 
to  be  real  if  j ^ i and  virtual  if  j = i . Virtual  transitions  are  represented  in  the  transi- 
tion matrix  by  nonzero  diagonal  elements.  This  formulation  allows  for  both  kinds  of 
transitions . 

Let  w^(t)dt  be  the  probability  of  leaving  state  i between  t and  t+dt  . The  quantity 
w^(t)  is  called  the  unconditional  waiting  time  density  and  is  expressed  in  terms  of  the 
hjj(.)  by: 

w (t)  = ^ h (t)  . (2) 

J ^ 

GENERATING  A SEMI-MARKOV  PROCESS 

There  are  two  general  types  of  random  mechanisms  that  give  rise  to  semi -Markov 
processes.  In  the  first  type  of  mechanism,  the  process  is  envisioned,  after  having  entered 
state  i , first  to  make  a random  draw  to  determine  the  successor  state  J , and  then  to 
make  another  random  draw  to  determine  the  amount  of  time  it  will  stay  in  state  i before 
going  to  state  J . The  outcome  of  the  first  draw  is  based  on  the  set  of  transition  probabil- 
ities p , 1 < i,  J < N . The  probability  p is  the  conditional  probability  that  the  next 


transition  is  to  state  J , given  that  the  last  transition  was  to  state  i . The  outcome  of 
the  second  draw  is  based  on  the  functions  • fhe  holding  time  density  functions  for 

duration  in  state  i given  that  the  next  transition  is  to  state  J.  Then: 


Ph 


and 


t,j(t)  = hy(t)/Py 


It  is  clear  from  the  definition  that  f,.(.)  integrates  to  unity;  from  equation  (1),  it  is 
seen  that:  ^ 


1 s i « N 


In  the  second  type  of  mechanism  the  process,  after  having  entered  state  i , makes 
N random  draws  -- one  draw  from  each  of  the  N holding  time  densities  • If 

then  determines  the  successor  state  and  length  of  time  in  state  i from  the  smallest  draw. 
Thus,  if  the  results  of  the  N draws  are  the  values  t^^  , 1 s j s N , and  if  t^^  = MjNftj^)  , 

then  the  next  transition  is  to  state  k and  the  length  of  time  the  process  holds  in  state  i 
before  going  to  k is  t^^  . 


An  example  of  the  first  type  of  mechanism  is  the  problem  of  the  whimsical  travelling 
salesman.  The  salesman  randomly  decides  which  city  he  will  visit  next,  given  that  he  is 
in  a particular  city  (decision  based  on  a random  draw  and  the  P^j's)  . Then  the  amount 

of  time  it  takes  the  salesman  to  travel  to  the  next  city  is  based  on  a random  draw  from  the 
density  of  travel  times  to  that  city. 


An  example  of  the  second  type  of  mechanism  is  the  reliability  of  a system.  If  a sys- 
tem is  composed  of  N critical  subsystems,  each  with  known  failure  time  densities,  then 
the  path  by  which  the  system  fails  and  the  time  to  failure  is  determined  by  the  smallest 
draw  from  the  failure  time  densities. 


by: 


The  holding  time  probabilities 


hy(.) 


for  the  two  types  of  mechanisms  are  expressed 


(3) 


(PiJ 

* f.. (1)11(1 -F  (t)) 

1 ^ m 


type  1 mechanism 
type  2 mechanism 


where 

0 

process  has  not  made 

Is  the  probability  that 
and  t+dt  . 


For  the  type  2 mechanism,  11  (1"F  . (t))  Is  the  probability  that  the 

kf^j 

a transition  to  a state  other  than  J by  time  t . Further,  fy(t)dt 
the  transition  to  j will  occur  In  the  next  Instant  of  time  between  t 


It  Is  highly  Important  In  treating  a semi -Markov  process  that  the  process  be  correctly 
Identified  as  either  type  1 or  2.  An  example  of  its  importance  is  when  there  is  a system 
composed  of  two  independent,  identical  subsystems  each  of  which  has  an  exponential  failure 

time  density  (l/r)e”*'^*^  . Then  the  mean  time  between  failure  is  r If  the  system  Is 
type  1,  but  1/2  r if  the  system  is  type  2. 


In  certain  cases,  a system  may  be  composed  of  some  states  that  obey  a type  1 transi- 
tion mechanism,  and  others  that  obey  a type  2 mechanism.  There  is  no  difficulty  in  model- 
ing such  a system,  providing  that  hy(t)  is  specified  properly  for  all  states. 

Reflection  shows  that  the  type  2 mechanism  is  common  and  is  capable  of  considerable 
generalization,  e.g.,  the  successor  state  nuiy  be  chosen  as  some  function  of  the  N draws, 
not  necessarily  the  state  yielding  the  smallest  draw.  The  possibility  of  anything  but  a type 
1 mechanism  is  seldom  mentioned  in  the  literature,  and  the  quantities  p^^  and  fy(.)  are 

usually  treated  as  fundamental,  with  li^j(*)  defined  in  terms  of  these.  Here,  h^j(.)  has 

been  taken  as  fundamental  to  show  that  it  can  be  expressed  in  various  ways  according  to 
the  way  in  which  the  semi -Markov  process  is  generated.  Other  mechanisms  for  choosing 
time  In  state  and  successor  state  can  be  found  and  give  rise  to  still  other  forms  for  h^^(.)  . 

Reference  4 discusses  the  case  in  which  the  choice  of  successor  state  is  probabilistic, 
conditional  on  the  value  of  the  waiting  time  in  the  current  sute. 


If  It  is  not  known  a priori  that  the  system  under  consideration  is  of  type  1,  type  2,  or 
otherwise,  it  may  be  possible  in  some  cases  to  construct  h^^(t)  directly  from  data.  When 

sufficient  data  is  available,  a convenient  r^resentatlon  of  h^^(t)  is  w^(t)  p|^(t)  , where 

w^(t)  is  the  waiting  time  density  for  transition  out  of  state  i , and  is  the  probability 
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that  the  system  will  make  the  next  transition  to  j , given  time  t and  current  state  i . 
This  latter  probability  may  be  easy  to  estimate  from  the  data. 

THE  SYSTEM  EQUATIONS 

This  subsection  derives  the  equations  obeyed  by  the  quantities  which  are  usually  of 
primary  interest  to  an  analyst.  No  attempt  at  completeness  is  made;  equations  for  other 
quantities  of  possible  interest  are  found  in  references  4,  5,  and  6. 


Let  0jj(t)  t>e  the  probability  that  the  system  is  in  state  j at  time  t , given  that  it 

entered  state  i at  time  zero.  The  quantities  are  the  state  probabilities  of  the 

semi -Markov  process.  Knowing  the  values  of  these  quantities  for  all  time  t is  equivalent 
to  knowing  the  probabilistic  time  behaviour  of  the  system.  The  governing  equation  for 
these  probabilities  is  now  derived.  A system  starting  in  state  i can  be  in  state  j at 
time  t in  the  following  ways: 

1.  i=j 


a.  System  never  leaves  state  i , or 

b.  System  leaves  state  1 but  returns  to  1 by  time  t . 

2.  ifi 

a . System  leaves  state  i and  manages  to  reach  j by  time  t . 


The  set  of  equations  that  describes  these  events  is: 

(^l^(t)  = + Z dT  h^(T)  (<j^j(t-T) 


(4) 


where 


t >0  , 1«  i,j  « N » , and  W^(t)  is  given  by: 

(t)  = 1 - /*  w (t)dt  . 

^ 0 ^ 


W 


(5) 


The  first  term  on  the  right-hand  side  of  equation  (4)  is  nonzero  when  i=j  , in  which  case 
it  represents  the  probability  that  the  system  did  not  leave  its  starting  state  i by  time  t . 
In  the  second  term,  the  quantity  hy^(T)dT  is  the  probability  that  the  system  will  make  a 


transition  to  state  k between  times  r and  t4<It 


This  is  multiplied  by  0.  .(t-r)  , the 


probability  that  the  system  will  proceed  to  state  J in  the  remaining  time  t-T  , having 
entered  state  k . The  product  is  Integrated  over  all  r between  0 and  t and  summed 
over  all  states  k . 
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When  N = 1 , equation  (4)  reduces  to  a single  (renewal)  equation;  for  N > 1 , it  is 
called  a "Markov  renewal  equation”  (references  5 and  6).  Many  quantities  of  interest  in 
the  analysis  of  semi-Markov  processes  obey  a Markov  renewal  equation;  this  Research 
Contribution  restricts  its  attention  to  such  quantities . 

When  the  semi-Markov  process  has  few  states  and  the  densities  are  simple  functions, 
it  is  sometimes  possible  to  take  the  Laplace  or  geometric  transform  of  equation  (4),  in- 
vert a matrix  and  take  the  inverse  transform  to  obtain  described  in  reference 

3.  Transform  techniques  are  usually  difficult  to  apply  and  a direct  numerical  evaluation 
of  equation  (4)  is  needed.  In  the  continuous -time  case,  there  is  little  literature  about  the 
numerical  solution  of  equation  (4).  Appendix  A derives  an  algorithm  for  solving  equation 
(4)  in  this  case.  In  the  discrete -time  case,  the  situation  is  somewhat  improved;  equation 
(4)  can  be  used  almost  "as  is"  to  produce  ^ recursive  manner.  This  is  also 

discussed  in  reference  4.  Appendix  B derives  an  algorithm  for  solving  equation  (4), 
which  has  excellent  numerical  accuracy  characteristics . 

Absorbing  States 


How  absorbing  states  are  treated  in  this  methodology  is  now  discussed.  An  absorbing 
state  is  a state  that  is  never  left  once  it  is  entered;  all  other  states  are  called  transient 
states . The  set  of  transient  states  is  denoted  T. 

In  brief,  the  effect  of  absorbing  states  is  to  simplify  the  solution  of  equation  (4). 

When  i is  an  absorbing  state,  it  is  seen  that: 

for  all  t and  all  j . Hence,  it  is  necessary  to  solve  equation  (4)  only  when  i is  a 
transient  state.  Also,  when  j is  also  a transient  state,  further  simplification  results 
from  noting  that  the  summation  in  equation  (4)  must  be  made  only  over  the  transient  states 
since  is  zero  whenever  k is  an  absorbing  state.  With  these  modifications,  the 

methodology  for  solving  equation  (4)  remains  effective  regardless  of  the  presence  or  ab- 
sence of  absorbing  states. 


Frequently,  the  analyst  is  particularly  interested  in  the  case  where  i is  a transient 
state  and  J is  an  absorbing  state.  In  this  case,  Is  s cumulative  distribution  func- 

tion, obtained  by  solving  equation  (4)  as  written.  However,  when  N is  large,  it  may  be 
computationally  advantageous  to  note  that  equation  (4)  can  be  rewritten  in  this  case  as: 


0,.(t)  « fdr  h.,(T)  £ j“dT  h,.  (T)  A ,(t-T) 

IJ  n y kcTO  ^ 


r 


$ 


(6) 


I 


I 

I 


where  the  summation  is  over  just  the  transient  states.  This  is  of  the  same  form  as 
equation  (4),  but  the  summation  is  over  fewer  states. 

Sometimes,  the  quantity  of  interest  is  the  density  function  corresponding  to 

the  distribution  function  d..(.)  . While  the  density  function  can  be  gotten  by  numerically 
differentiating  d. .(•)  f ft  vJill  usually  be  better  practice  to  differentiate  both  sides  of 
equation  (6)  and  obtain  the  following  equation  for  : 

d*..(t)  = h (t)  + L f ClT  h (t)  d*.  .(t-T)  . (7) 

1]  IJ  Q IK  KJ 


First  Passage  Time  Distributions 


Sometimes,  the  quantity  of  main  interest  is  the  distribution  of  the  "first  passage 
time,”  i.e.,  the  time  to  the  first  visit  of  state  j given  start  in  state  i . A straight- 
forward way  to  calculate  this  distribution  is  to  make  state  j absorbing,  then  apply  the 
above  methodology  to  calculate  0y(*)  or  d*j.(.)  • 


Expected  Time  in  State 

Sometimes  the  quantity  of  interest  is  the  expected  amount  of  time  spent  in  a given 
state  in  a given  time  interval  [O,  T]  , The  random  variable  X..  is  defined  equal  to  one 

if  the  process  is  in  state  j at  time  t , given  start  in  state  i at  time  zero,  and  equal  to 
zero  otherwise.  Then  the  total  time  spent  in  state  j is: 

jV(t)dt 

0 J 

and  the  expected  time  in  state  j is: 


rx  (t)]dt  = /'^d„(t)dt  . (8) 

0 ^ 0 ^ 


Hence,  the  desired  expectation  is  gotten  by  solving  equation  (4)  for  d.  .(.)  * then  inte 
grating  it. 


A CONTINUOUS -TIME  APPLICATION 


A simple  stochastic  process  will  now  be  discussed.  Also,  the  equations  developed 
in  the  preceding  section  will  be  applied  to  produce  the  numerical  solution  of  the  probabil- 
istic time  behaviour  of  the  process.  The  numerical  solution  will  then  be  compared  with 
the  analytic  solution. 

The  problem  to  be  considered  is  of  an  aircraft  searching  for  a submerged  submarine 
in  the  ocean.  The  aircraft  drops  sensors  to  detect  the  submarine.  The  aircraft  is  as- 
sumed to  detect  the  submarine  at  a constant  rate  X2  , where  reciprocal  of  the 

mean  time  to  detect.  After  the  aircraft  detects  the  submarine,  it  drops  additional  sensors 
to  localize  (refine  the  location)  the  submarine.  Localization  is  assumed  to  take  place  at 
a constant  rate  X^  , where  X^  is  then  reciprocal  of  the  mean  time  to  localize.  Once 

the  aircraft  has  localized  the  submarine,  the  mission  is  considered  successful  and  the 
aircraft  returns  to  base. 

While  the  aircraft  is  attempting  localization,  however,  the  submarint  may  move  out 
of  detection  range  of  the  sensors  and  the  aircraft  may  lose  contact.  The  aircraft  is  as- 
sumed to  lose  contact  on  the  submarine  while  attempting  to  localize  it  at  a constant  rate 
\ , where  X is  the  reciprocal  of  the  mean  time  to  lose  contact . If  the  aircraft  loses 

contact  on  the  submarine,  additional  sensors  are  laid  to  redetect  the  submarine.  The 
redetectlon  rate  is  assumed  to  be  equal  to  the  detection  rate  X^  . Given  that  the  aircraft 

can  search  for  a long  time  and  has  an  inexhaustible  supply  of  sensors,  then  the  probabili- 
tiej  of  the  aircraft  being  in  the  search,  detection,  and  localization  states  as  a function  of 
time  must  be  determined. 

The  elementary  process  just  described  is  a three -state  system  consisting  of  two 
transient  states  and  one  absorbing  state  (figure  1).  The  arrows  in  the  figure  indicate 
the  direction  of  transitions  from  each  of  the  system. 


FIG.  1:  THREE-STATE  STOCHASTIC  PROCESS 
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State  1 (localization)  is  the  absorbing  state  because  once  entered  it  is  never  left,  while 
states  2 and  3 (detection  and  search,  respectively)  are  transient  states.  The  holding  time 
densities  between  states  of  the  system  are  exponential  and  are  indicated  in  the  figure. 


The  system  in  figure  1 is  assumed  to  make  transitions  according  to  a type  2 mechan- 
ism. The  localization  state  and  the  search  state  compete  for  transitions  from  the  detection 
state  on  a time  basis.  Using  equation  (3),  the  probability  transition  matrix  is: 


1 

1 /6(t— ) 


2 

0 


H(t)  = 2|  0 


0 


1 

V 2 


(9) 


where  6(t-<»)  is  a formalism  to  mean  that  state  1 is  never  left,  once  entered  (except  at 
t=«)  . From  equation  (5),  the  W(t)  matrix  is  given  by: 


W(t)=|  e’^^l''’^3^'' 


2 


(10) 


Using  equations  (9)  and  (10)  in  the  recursive  equation  (A -5)  of  appendix  A,  •(!)  is 
produced  numerically.  Figures  2 and  3 present  the  solution  ?(t)  for  \j=l/24  , 

X2=1/12  , X2=1/18  , and  At=0.5  . Note  that  * *i2~*13^  for  all  time  because 

state  1 is  the  absorbing  state.  Figure  2 displays  the  solution  if  the  system  initially 
started  in  state  3.  The  probability  of  being  in  the  absorbing  state  increases  rapidly 

and  then  asymptotes  to  1.  The  probability  of  being  in  state  3,  given  that  the  process 
started  there  f . , rapidly  decays  to  zero.  The  probability  of  being  in  state  2 gflven 

that  the  system  started  in  state  3 » buiW®  up  0 at  time  zero  to  a maximum  of 

.39  at  time  14  and  then  decays  to  zero.  Figure  2 shows  similar  behaviour  when  the 
system  starts  in  state  2. 

Table  1 compares  the  numerical  solution  of  equation  (A -5)  for  and  # the 

absorption  probabilities,  with  a numerical  evaluation  of  the  analytic  solution  of  the  system 
of  differential  equations: 
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FIG.  3:  EXAMPLE  SOLUTION  FOR  SYSTEM  STARTING  IN  STATE  2 


i 


11 


d? 

IT  " '2*3  ^3*2 


^=-(VV*2-^S’3 


d$ 

— =-  = > $ 
dt  12 


(11) 


which  also  describes  the  process  in  figure  1.  Here  f , i , and  $ represent  the 

X ^ <3 

probabilities  of  being  in  states  1,  2,  and  3,  respectively.  The  general  solution  of  equation 
(11)  yields  equations  with  undetermined  constants.  These  constants  are  evaluating  by  im- 
posing initial  conditions,  i.e.,  4-(0)=l,  $„(0)=5  (0)=0  for  starting  in  state  3,  and  ? (0)=1  , 

o ^ L 2 

$2(0)=$j^(0)=0  for  starting  in  state  2. 


TABLE  1 

COMPARISON  OF  NUMERICAL  SOLUTION  WITH  EXACT  SOLUTION 
Numerical  solution  Exact  solution® 

Time  (units)  *31 *21  *31  *21 


0 

0 

0 

0 

0 

5 

. 03264 

.16792 

.03261 

. 16786 

10 

. 10081 

. 28325 

.10074 

. 28315 

15 

.17953 

.37072 

.17941 

.37059 

20 

.25805 

.44194 

.25788 

.44176 

25 

. 33203 

.50251 

.33182 

. 50231 

30 

. 39998 

. 55534 

.39973 

. 55511 

35 

.46162 

.60203 

.46133 

.60177 

40 

.51721 

. 64358 

. 51688 

.64329 

45 

. 56719 

.68069 

. 56683 

.68037 

50 

.61206 

.71389 

.61167 

.71355 

100 

.87058 

.90462 

.86996 

.90411 

a 


Numerical  evaluation  of  analytic  solution. 


As  can  be  seen  in  table  1 the  difference  between  the  numerical  solution  and  the  exact 
solution  is  less  than  1 part  of  10^.  A smaller  step  size  would  have  yielded  better  accur- 
acy but  required  much  more  computation  time.  The  numerical  error  inherent  in  the 
trapezoidal  integration  rule  is  given  by  equation  (A -3)  as: 
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A DISCRETE -TIME  APPLICATION 

This  section  models  a dogfight  between  "Snoopy"  and  the  "Red  Baron"  as  a semi-Markov 
process  in  discrete  time  to  obtain  the  probability  that  Snoopy  wins  as  a function  of  time. 

The  numerical  solution  is  then  compared  with  the  analytic  solution. 

The  dogfight  is  thought  of  as  being  in  one  of  five  possible  states  at  any  time;  (1)  win 
for  Snoopy,  (2)  offensive  position  for  Snoopy,  (3)  neutral,  (4)  defensive  position  for  Snoopy, 
and  (5)  loss  for  Snoopy.  The  transitions  between  states  are  governed  by  the  transition 
matrix  below; 


Modeling  the  dogfight  as  a type  1 mechanism  seems  most  sensible,  and  the  holding  time 
densities  are  defined  by; 

f (n)  = p q/'  \ n=l,2, ...  fori=2,3,4,  withp  =p  =.2,  p_=.l,  q.=l-p.  . 

Ijli 

States  1 and  5 are  absorbing  states  so  that  their  holding  time  densities  need  not  be  speci- 
fied. This  choice  of  fy(«)  fs  made  primarily  for  analytic  convenience  and  these  functions 

do  not  represent  realistic  holding  time  distributions  in  actual  dogfights . 

Of  interest  here  is  the  probability  that  Snoopy  wins,  as  a function  of  time.  To  make 
matters  concrete,  it  is  assumed  that  the  Red  Baron  ambushes  Snoopy  (the  engagement 
starts  in  state  4),  and  also  that  the  fight  breaks  off  (due  to  fuel  exhaustion)  after  2 minutes 
if  no  one  wins  before  that  time . 

The  discrete-time  version  of  equation  (7)  is  wanted;  specifically,  desired 

since  state  1 is  the  absorbing  state  corresponding  to  a win  by  Snoopy,  and  the  fight  starts 

in  state  4.  Also,  At  » 1 second,  so  that  t »n,  n^O, ...,  120  . Then; 

n 

n 4 

0*„(n)  ■ h (n)  + E E h k(T)0*  (n-T)  (12) 

T-1  k«2 


I 
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just  as  equation  (B-1)  in  appendix  B was  gotten  from  equation  (4).  To  get  , it  is 

necessary  to  solve  for  all  0*^  , 1=2, 3,4  . To  get  the  final  form  of  the  recursion 
scheme,  equation  (12)  is  rewritten  in  matrix/vector  form  as; 

n 

0*(n)  = h(n)+  G(T)0*(n-T)  (13) 

T=1 


where  j0*(n)  is  the  3x1  vector  with  components  » 1=2,3, 4;  h(n)  is  the  corres- 

ponding vector  with  components  h^j(n)  , and  G is  the  3x3  submatrix  of  H , consisting 

only  of  the  transient  states  2,  3,  and  4.  Equation  (13)  is  a discrete  form  of  equation  (7). 
Alternatively,  0*^^  could  be  gotten  by  solving  equation  (B-1)  for  0^^  and  "differentiating” 

to  get  0*^^  . But  this  way  is  computationally  superior  in  the  same  way  that  equation  (7) 

is  superior  to  equation  (6)  when  absorbing  states  are  involved. 

Equation  (13)  can  be  solved  using  the  z -transform  (reference  4),  and  the  probability 
of  a win  by  Snoopy  given  the  start  in  state  4 is: 


(•'®)  (*'4) 


1 ^vH-l 

- .09(,8) 


n 2 1 . 


A comparison  of  the  exact  values  and  the  values  gotten  by  recursive  solution  using  equa- 
tion (13)  is  given  in  table  2. 


There  Is  agreement,  in  the  table,  through  the  first  nine  digits.  However,  it  is  possible 
that  the  numerical  solution  is  more  accurate  than  the  "exact”  solution,  depending  on  how 
well  the  CNA  computer  raises  numbers  to  large  powers. 

Figure  4 gives  the  cumulative  probability  of  a win  for  Snoopy  as  a function  of  the  length 
of  the  dogfight  for  all  three  surting  states,  2,  3,  and  4.  The  probability  of  a loss  for 
Snoopy  could  be  calculated  similarly  by  changing  the  vector  h(n)  to  the  vector  with  com- 
ponents • The  exchange  ratio  as  a function  of  time  can  be  gotten  by  dividing  the 

probability  M a win  by  the  probability  of  a loss . 
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Probability  of  win 


TABLE  2 


COMPARISON  OF  EXACT  VALUES  AND  VALUES 
FROM  RECURSIVE  SOLUTION 


Time 

Numerical 

Exact 

(sec) 

solution 

solution 

0 

0 

0 

10 

9.274338537F“ 

3 

9. 27433054;?"3 

20 

1,021862002;?" 

2 

1.021862002;?" 

30 

7. 384697292F" 

3 

7.384697295;?" 

40 

4,980723436K" 

3 

4,980723439;?" 

50 

3. 319363533K" 

3 

3.319363535;?" 

60 

2.207700531.*?" 

3 

2.207700533;?" 

70 

1.467846696A" 

3 

1.467846698;?" 

60 

9.758826697P" 

4 

9. 75882671K"4 

90 

6. 487997539K" 

4 

6.487997548;?" 

100 

4.313433694K’" 

4 

4, 313433701^" 

110 

2.867711621;?" 

4 

2.867711626”" 

120 

1.90654829P“4 

1.90654R293R" 

a 


Numerical  evaluation  of  analytic  solution. 


FIQ,  4:  PROBABILITY  OF  WIN  FOR  SNOOPY  vt.  DOGFIGHT  DURATION 
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DISCUSSION  OF  ALGORITHMS 


In  the  discrete -time  case,  algorithms  for  the  solution  of  equation  (4)  are  easy  to 
write.  The  algorithm  in  appendix  B is  fast  running  and  suitable  for  time -sharing  usage. 
The  answers  are  exact  except  for  truncation  and  roundoff  errors,  both  of  which  accumu- 
late slowly. 

In  the  continuous -time  case,  the  situation  is  not  as  good.  An  algorithm  based  on  the 
trapezoidal  rule  is  given  in  appendix  A.  It  is  easy  to  program,  but  requires  a small  step 
size  for  accuracy.  The  running  time  for  the  algorithm  is  proportional  to  n^  , where  n 
is  the  number  of  subdivisions  of  the  time  axis,  so  that  small  step  sizes  imply  long  running 
times.  Appendix  A also  contains  an  algorithm  based  on  Simpson's  rule.  This  algorithm 
gives  about  three  more  digits  of  accuracy  than  the  trapezoidal  algorithm  for  the  same  step 
size  for  most  applications;  however,  the  Simpson  algorithm  is  more  complicated  than  the 
trapezoidal  algorithm  and  is  inherently  slower  running.  All  of  the  above  algorithms  can 
be  programmed  in  various  ways,  e.g.,  to  emphasize  speed  or  to  minimize  storage. 

Relatively  little  literature  exists  concerning  the  numerical  solution  to  equation  (4). 
Many  authors  suggest  the  use  of  transform  methods,  involving  a numerical  inversion  of 
the  transform.  This  is  a nontrivial  operation  in  most  cases.  The  idea  of  a direct  numer- 
ical attack  on  equation  (4)  is  not  new;  reference  7,  for  example,  outlines  an  algorithmic 
procedure  very  similar  to  that  used  in  appendix  A.  although  no  discussion  of  numerical 
accuracy  or  running  time  is  given.  Reference  8 describes  a solution  based  on  expansion 
in  power  series,  resulting  in  a simple  recursive  scheme  for  determining  the  unknown 
coefficients.  The  technique  must  be  tailored  to  the  specific  distributions  in  use,  however. 
Reference  9 discusses  a rather  elaborate  scheme  using  spline  functions  and  Lobatto  inte- 
gration to  do  calculations  similar  to  those  needed  for  equation  (4) . But  it  is  not  clear  that 
such  a complex  technique  is  needed  for  routine  solution  of  equation  (4). 

Reference  10  gives  an  annotated  bibliography  on  the  computational  aspects  of  proba- 
bilistic modeling.  In  it  will  be  found  very  few  references  to  the  computational  problems 
involved  in  solving  Markov  renewal  equations  (equation  (4)).  It  seems  clear  that  consider- 
ably more  algorithmic  research  is  needed  before  such  equations  can  be  routinely  solved 
in  a fast  and  accurate  manner. 
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APPENDIX  A 

CONTINUOUS -TIME  SOLUTION  FOR  THE  SYSTEM  EQUATIONS 

In  this  appendix,  two  algorithms  are  derived  for  solving  equation  (4)  (from  main 
text)  by  numerical  integration.  First,  an  algorithm  is  derived  based  on  the  trapezoidal 
rule;  then  an  algorithm  is  derived  based  on  Simpson’s  rule. 

TRAPEZOIDAL  RULE 

Equation  (4)  is  solved  by  finding  on  a set  of  equally  spaced  points  in  time. 

The  time  points  are  denoted  size,  i.e,,  the  interval  between 

the  time  points,  is  denoted  At  . The  time  points  can  be  represented  as; 

t = At*n,  0 s n s m . 
n 

With  0..(O)  = 6 , the  solution  for  t , n >0  , is  now  desired, 

ij  IJ  n 

For  convenience,  equation  (4)  is  rewritten  in  matrix  notation  as: 


0(t  ) = W(t_,)  + J " dT  H(t)  0(t^-T) 

n n n 


(A-l) 


where 


0(V  = 

H(t„)  = 


1 s i,J  s N 


and  W(t  ) is  the  diagonal  matrix  with  elements  W.(t  ) down  the  diagonal  and  zeros 
n in 

elsewhere.  It  is  observed  that  when  t=0  , the  argument  of  ^ in  the  integrand  is  the 

same  as  that  of  0 on  the  left  of  the  equation.  To  capitalize  on  this,  equation  (A-l)  is 

rewritten  as  follows: 


0(t  ) = W(t  ) + J " dr  H(t)  0(t  -t  ) 

" " 0 " 

t 

+ J*  " dr  H(t)  0(t  -t)  . 


(A -2) 


A-l 


Equation  (A -2)  is  evaluated  using  the  trapezoidal  integration  rule  (reference  11),  i.e., 


J i(y)  dy  = ^ Cf(x^)  + f(x2)l 


Ax 


(A -3) 


where  x^  and  x^  are  separated  by  an  interval  Ax  , f(.)  is  the  function  to  be  inte- 
grated, and  the  last  term  on  the  right  is  the  error  term.  Applying  equation  (A -3)  without 
the  error  term  produces; 

0(tn)  = W(t^)  + ^ rH(O)0(t^)  + H(t^)0(t^-t^)] 

n 

+ At  p H(y0(t^-tj^)  (A -4) 

[H(t,)0(t  -t,)  + H(t  )0(O)1 
z 1 n 1 n 

where  H(0)  = H(tQ)  has  been  used.  Solving  for  produces; 

0(t^)  = [I  H(0)l"^  X rW(y  + At  X]  "T  * <^-5) 

k=l 

This  equation  is  the  basic  recursive  scheme  used  numerically  to  produce  the  solution  of 
equation  (4)  on  the  points  t^,  1 < n < m . The  solution  is  started  with  0(0)  = W(0)  = 1 , 

and  the  solution  for  0(1^)  is  gotten  by  applying  equation  (A -5),  etc.  In  this  way,  the 

solution  of  equation  (4)  for  any  finite  t^  is  produced. 

The  matrix  [I  “ H (0)  1 In  equation  (A  -5)  need  be  inverted  only  once  at  the  start 

of  the  numerical  scheme.  If  the  transition  matrix  H(t)  contains  only  density  functions 
that  are  zero  at  time  zero,  no  inversion  is  required. 

The  transition  matrix  should  contain  only  density  functions  that  are  continuous  in  the 
time  interval  of  iitterest.  Discontinuous  density  functions  may  cause  large  numerical 
Inaccuracies  in  the  numerical  scheme. 


A -2 


At  ea;:h  step,  the  solution  of  equation  (A -5)  must  satisfy  the  following  condition: 


N 

E 0 (t  ) = 1 . for  1 « i =£  N 

4=1  ” 


(A -6) 


This  condition  serves  to  check  on  both  the  numerical  accuracy  and  stability  of  the  solution. 

The  relative  numerical  accuracy  of  the  iterative  scheme  is  one  order  of  magnitude 
less  than  the  numerical  accuracy  of  the  trapezoidal  integration  rule.  This  occurs  because 
of  the  compounding  of  numerical  errors  in  the  iterative  scheme.  The  relative  error  for 
the  trapezoidal  rule  is  given  by: 

_(Ax)^  r(g) 

12  f(?) 

where  f(.)  is  the  function  being  integrated. 


SIMPSON'S  RULE 


In  this  subsection,  Simpson  rules  are  applied  to  equation  (A-1).  Simpson's  rule  and 
Simpson's  3/8  rule  (reference  11)  are  given  by: 


f(t)dt  = ^ (fg  + ^1+^2^  ■ 


'0 


and 


90 


_ lAL  ...  ^ - 3(At)^ 


; f(t)dt  = ^ (fQ+afj+Sf^+fg)  - 


80 


(A-7) 

(A -8) 


respectively,  where  At  is  the  step  size  and  f^  = f(tj)  . For  larger  intervals  of  integra- 
tion, equation  (A-7)  is  applied  repeatedly  to  obtain: 


‘o 


(A^) 


The  right-hand  side  is  denoted  by  S(tg,t^)  . Equation  (A -9)  is  also  known  as  Simpson's 
rule.  Note  that  n must  be  even  and  greater  than  or  equal  to  two. 


These  rules  are  applied  for  the  numerical  solution  of  equation  (A-1),  and  it  will  be 
seen  that  the  cases  "n  even”  and  ”n  odd”  must  be  handled  differently. 

A -3 


Case  I:  n Even,  n s 2 


The  integral  in  equation  (A-1)  is  split  at  the  point  t2  , obtaining; 

0(t  ) = W(t„)  + ^ [H(O)0(t^)  + 4H(t,)0(t  ,)  + H(t  J0(t^  J]  + S(t,.  t ) 

n n o n l n-i  z n-l  i n 

or  (A-10) 

0(y  = [I  - H(0)]"^  fW(y  + ^ C4H(tj),rf(t^_^)  + H(t2))i(t^_2)]  + t^)}  . 


Case  II;  n Odd,  n i 3 


The  integral  in  equation  (A-1)  is  split  at  t^  , obtaining; 


3At 


d(t^)  = W(t^)  + ^ rH(O)0(t^)  + 3H(t^)0(t^_^)  + 3H(t2)0(t^.2)  + H(t3)0(t^_3)l 


+ S(t  ,t  ) 
j n 


or 


0(t^)  = [I  - ^ H(0)]"^  nV(t^)  + ^ [3H(tj)0(t^_j)  + 3H(t2)0(t^.2> 


(A-11) 


Equations  (A-10)  and  (A-11)  are  the  counterpart  of  equation  (A -5).  It  should  be  noted  that 
there  is  no  way  to  get  0(t  ) by  these  formulas.  This  value  must  be  input  independently, 
if  the  analyst  wishes  to  get  the  increased  accuracy  that  Simpson’s  rule  provides  over  the 
trapezoidal  rule.  Indeed,  as  the  analyst  passes  to  more  accurate  numerical  integration 
schemes,  he  will  have  to  input  increasingly  more  points  to  begin  the  calculation.  In 
practice,  this  creates  no  problem  --  the  analyst  might  use  something  like  the  trapezoidal 
rule  to  get  these  initial  points,  then  shift  to  the  more  accurate  scheme. 
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DISCRETE -TIME  SOLUTION  FOR  THE  SYSTEM  EQUATIONS 

When  the  holding  time  variables  are  discrete  and  take  their  values  on  the  lattice 
points  t^  = 4fn,  Osnsm  ,asin  appendix  A,  then  equation  (4)  can  be  solved  with  a 

simple  recursive  algorithm  which  has  excellent  numerical  error  characteristics.  The 
algorithm  is  derived  in  reference  4,  but  is  rederived  here  for  completeness. 

In  the  discrete  case,  equation  (A-1)  in  appendix  A reduces  to: 


t 

n 


0(tn)  = W(t  ) + Yi  H(t)  0(t  -t) 
n n n 


However,  it  has  been  required  that  the  holding  time  densities  h.j(.)  have  no  impulse 
component  at  the  origin,  so  that  in  the  discrete  case  H(0)  is  the  zero  matrix.  Then: 


n 

0(t  ) = W(t  )+  i;H(T)0(t  -t) 
n n n 


(B-1) 


f It  is  also  seen  that  ^(t  ) has  been  ejqiressed  in  terms  of  <j>(t, ),  0 ^k^n-l.  Starting  with 

I n K 

1^(1^)  = W(tQ),  equation  (B-1)  permits  the  recursive  computation  of  (t^)  . 

Equation  (B-1)  is  remarkably  well  suited  to  numerical  computation  since  (1)  there 
are  no  matrix  inversions,  (2)  there  are  no  subtractions  or  divisions,  and  (3)  there  are 
no  scaling  problems  since  all  the  quantities  involved  are  probabilities  and  hence  all 
i scaled  between  zero  and  one.  It  may  therefore  be  anticipated  that  roundoff  and  trunca- 

tion errors  will  accumulate  very  slowly,  even  in  problems  involving  large  dimensions, 
i 
i 


